Load the necessary libraries
library(rstanarm) #for fitting models in STAN
library(brms) #for fitting models in STAN
library(coda) #for diagnostics
library(bayesplot) #for diagnostics
library(ggmcmc) #for MCMC diagnostics
library(rstan) #for interfacing with STAN
library(emmeans) #for marginal means etc
library(broom) #for tidying outputs
library(DHARMa) #for residual diagnostics
library(tidybayes) #for more tidying outputs
library(ggeffects) #for partial plots
library(broom.mixed)#for tidying MCMC outputs
library(tidyverse) #for data wrangling etc
library(patchwork) #for multiple plots
library(ggridges) #for ridge plots
source('helperFunctions.R')
To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.
A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.
Barramundi
Format of norin.csv data files
| FISHID | MASS | TRIAL | SMR_contr | CHANGE |
|---|---|---|---|---|
| 1 | 35.69 | LowSalinity | 5.85 | -31.92 |
| 2 | 33.84 | LowSalinity | 6.53 | 2.52 |
| 3 | 37.78 | LowSalinity | 5.66 | -6.28 |
| .. | .. | .. | .. | .. |
| 1 | 36.80 | HighTemperature | 5.85 | 18.32 |
| 2 | 34.98 | HighTemperature | 6.53 | 19.06 |
| 3 | 38.38 | HighTemperature | 5.66 | 19.03 |
| .. | .. | .. | .. | .. |
| 1 | 45.06 | Hypoxia | 5.85 | -18.61 |
| 2 | 43.51 | Hypoxia | 6.53 | -5.37 |
| 3 | 45.11 | Hypoxia | 5.66 | -13.95 |
| FISHID | Categorical listing of the individual fish that are repeatedly sampled |
| MASS | Mass (g) of barramundi. Covariate in analysis |
| TRIAL | Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen. |
| SMR_contr | Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia) |
| CHANGE | Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'. |
norin <- read_csv('../public/data/norin.csv', trim_ws=TRUE)
## Rows: 180 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): TRIAL
## dbl (4): FISHID, MASS, SMR_contr, CHANGE
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(norin)
## Rows: 180
## Columns: 5
## $ FISHID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ MASS <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.9…
## $ TRIAL <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity", …
## $ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589, …
## $ CHANGE <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.0…
We should start by declaring any categorical variables (including an varying effects) as predictors.
Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i}\\ \beta_0 \sim{} \mathcal{N}(16, 35)\\ \beta_{1-6} \sim{} \mathcal{N}(0, 70)\\ \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of temperature and (centered) mean fish size on SDA peak. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual fish.
Lets begin by declaring the categorical predictors and random effects as factors.
norin <- norin %>% mutate(FISHID=factor(FISHID),
TRIAL=factor(TRIAL))
Each of the fish were exposed to each of the three trials (treatments). So lets explore the distributions of responses within each of these trials.
ggplot(norin, aes(y=CHANGE, x=TRIAL)) + geom_boxplot()
Conclusions:
The researchers considered that physiological plasticity might be effected by the fishes basal metabolic rate. For example, a fish with relatively high metabolism might have less scope to increase this metabolism than a fish with a lower metabolism. Therefore, the SMR under control conditions (just prior to the specific trial) could be added as a continuous covariate.
Doing so would introduce two additional assumptions:
ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
geom_smooth(method='lm') + geom_point()
## `geom_smooth()` using formula 'y ~ x'
ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
geom_smooth() + geom_point()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
It might also be worth exploring how consistent the trial effect is across the different fish. This can give us an idea of whether the addition of a random intercept/slope model would be useful.
ggplot(norin, aes(y=CHANGE, x=as.numeric(FISHID), color=TRIAL)) +
geom_point() + geom_line()
Conclusions:
Finally, as an acknowledgement that the metabolic response might be influenced
by the mass of the fish, the researchers contemplated including fish Mass as a continuous covariate. There are numerous alternative ways to incorporate covariates that are known to impact on a response.
Lets explore the relationship between the response and Mass separately for each Trial.
ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) +
geom_point() +
geom_smooth(method='lm')
## `geom_smooth()` using formula 'y ~ x'
Conclusions:
In rstanarm, the default priors are designed to be
weakly informative. They are chosen to provide moderate regularlization
(to help prevent overfitting) and help stabalise the computations.
norin.rstanarm <- stan_glmer(CHANGE ~ (1|FISHID)+TRIAL*SMR_contr+MASS,
data = norin,
family = gaussian(),
iter = 5000,
warmup = 2000,
chains = 3,
thin = 5,
refresh = 0)
norin.rstanarm %>% prior_summary()
## Priors for model '.'
## ------
## Intercept (after predictors centered)
## Specified prior:
## ~ normal(location = 20, scale = 2.5)
## Adjusted prior:
## ~ normal(location = 20, scale = 85)
##
## Coefficients
## Specified prior:
## ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
## Adjusted prior:
## ~ normal(location = [0,0,0,...], scale = [178.99,178.99,135.15,...])
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.03)
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details
This tells us:
mean(norin$CHANGE)) and
the scaled standard deviation of 85 is based on multiplying the scaling
factor by the standard deviation of the response.2.5*sd(norin$CHANGE)
## [1] 84.6135
2.5*sd(norin$CHANGE)/apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
## TRIALHypoxia TRIALLowSalinity
## 178.99305 178.99305
## SMR_contr MASS
## 135.14517 10.59372
## TRIALHypoxia:SMR_contr TRIALLowSalinity:SMR_contr
## 34.19156 34.19156
rstanarm, this is a exponential with a rate
of 1 which is then adjusted by devision with the standard deviation of
the response.1/sd(norin$CHANGE)
## [1] 0.02954611
Lets now run with priors only so that we can explore the range of values they allow in the posteriors.
norin.rstanarm1 <- update(norin.rstanarm, prior_PD=TRUE)
norin.rstanarm1 %>%
ggpredict(~SMR_contr*TRIAL) %>%
plot(add.data=TRUE)
## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
#OR
norin.rstanarm1 %>%
ggemmeans(~SMR_contr*TRIAL) %>%
plot(add.data=TRUE)
Conclusions:
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):
median(norin$CHANGE)mad(norin$CHANGE)sd(norin$CHANGE) / apply(model.matrix(~scale(SMR_contr)*TRIAL + scale(MASS), norin), 2, sd)1 / sd(norin$CHANGE)I will also overlay the raw data for comparison.
norin.rstanarm2 <- stan_glmer(CHANGE ~ (1|FISHID)+TRIAL*scale(SMR_contr)+scale(MASS),
data = norin,
family = gaussian(),
prior_intercept = normal(17, 35, autoscale = FALSE),
prior = normal(0, 70, autoscale = FALSE),
prior_aux=rstanarm::exponential(0.03, autoscale = FALSE),
prior_covariance = decov(1, 1, 1, 1),
prior_PD = TRUE,
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)
norin.rstanarm2 %>%
ggpredict(~SMR_contr * TRIAL) %>%
plot(add.data = TRUE)
## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
Now lets refit, conditioning on the data.
norin.rstanarm3 <- update(norin.rstanarm2, prior_PD=FALSE)
get_variables(norin.rstanarm3)
## [1] "(Intercept)"
## [2] "TRIALHypoxia"
## [3] "TRIALLowSalinity"
## [4] "scale(SMR_contr)"
## [5] "scale(MASS)"
## [6] "TRIALHypoxia:scale(SMR_contr)"
## [7] "TRIALLowSalinity:scale(SMR_contr)"
## [8] "b[(Intercept) FISHID:1]"
## [9] "b[(Intercept) FISHID:2]"
## [10] "b[(Intercept) FISHID:3]"
## [11] "b[(Intercept) FISHID:4]"
## [12] "b[(Intercept) FISHID:5]"
## [13] "b[(Intercept) FISHID:6]"
## [14] "b[(Intercept) FISHID:7]"
## [15] "b[(Intercept) FISHID:8]"
## [16] "b[(Intercept) FISHID:9]"
## [17] "b[(Intercept) FISHID:10]"
## [18] "b[(Intercept) FISHID:11]"
## [19] "b[(Intercept) FISHID:12]"
## [20] "b[(Intercept) FISHID:13]"
## [21] "b[(Intercept) FISHID:14]"
## [22] "b[(Intercept) FISHID:15]"
## [23] "b[(Intercept) FISHID:16]"
## [24] "b[(Intercept) FISHID:17]"
## [25] "b[(Intercept) FISHID:18]"
## [26] "b[(Intercept) FISHID:19]"
## [27] "b[(Intercept) FISHID:20]"
## [28] "b[(Intercept) FISHID:21]"
## [29] "b[(Intercept) FISHID:22]"
## [30] "b[(Intercept) FISHID:23]"
## [31] "b[(Intercept) FISHID:24]"
## [32] "b[(Intercept) FISHID:25]"
## [33] "b[(Intercept) FISHID:26]"
## [34] "b[(Intercept) FISHID:27]"
## [35] "b[(Intercept) FISHID:28]"
## [36] "b[(Intercept) FISHID:29]"
## [37] "b[(Intercept) FISHID:30]"
## [38] "b[(Intercept) FISHID:31]"
## [39] "b[(Intercept) FISHID:32]"
## [40] "b[(Intercept) FISHID:33]"
## [41] "b[(Intercept) FISHID:34]"
## [42] "b[(Intercept) FISHID:35]"
## [43] "b[(Intercept) FISHID:36]"
## [44] "b[(Intercept) FISHID:37]"
## [45] "b[(Intercept) FISHID:38]"
## [46] "b[(Intercept) FISHID:39]"
## [47] "b[(Intercept) FISHID:40]"
## [48] "b[(Intercept) FISHID:41]"
## [49] "b[(Intercept) FISHID:42]"
## [50] "b[(Intercept) FISHID:43]"
## [51] "b[(Intercept) FISHID:44]"
## [52] "b[(Intercept) FISHID:45]"
## [53] "b[(Intercept) FISHID:46]"
## [54] "b[(Intercept) FISHID:47]"
## [55] "b[(Intercept) FISHID:48]"
## [56] "b[(Intercept) FISHID:49]"
## [57] "b[(Intercept) FISHID:50]"
## [58] "b[(Intercept) FISHID:51]"
## [59] "b[(Intercept) FISHID:52]"
## [60] "b[(Intercept) FISHID:53]"
## [61] "b[(Intercept) FISHID:54]"
## [62] "b[(Intercept) FISHID:55]"
## [63] "b[(Intercept) FISHID:56]"
## [64] "b[(Intercept) FISHID:57]"
## [65] "b[(Intercept) FISHID:58]"
## [66] "b[(Intercept) FISHID:59]"
## [67] "b[(Intercept) FISHID:60]"
## [68] "sigma"
## [69] "Sigma[FISHID:(Intercept),(Intercept)]"
## [70] "accept_stat__"
## [71] "stepsize__"
## [72] "treedepth__"
## [73] "n_leapfrog__"
## [74] "divergent__"
## [75] "energy__"
posterior_vs_prior(norin.rstanarm3, color_by='vs', group_by=TRUE,
facet_args=list(scales='free_y'),
regex_pars = "^.Intercept|TRIAL|SMR_contr|MASS|sigma|Sigma")
##
## Drawing from prior...
Conclusions:
norin.rstanarm3 %>%
ggpredict(~SMR_contr * TRIAL) %>%
plot(add.data = TRUE)
## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
##OR
norin.rstanarm3 %>%
ggemmeans(~SMR_contr * TRIAL) %>%
plot(add.data = TRUE)
In brms, the default priors are designed to be weakly
informative. They are chosen to provide moderate regularlization (to
help prevent overfitting) and help stabalise the computations.
Unlike rstanarm, brms models must be
compiled before they start sampling. For most models, the compilation of
the stan code takes around 45 seconds.
norin.form <- bf(CHANGE ~ (1|FISHID)+TRIAL*SMR_contr+MASS,
family = gaussian()
)
options(width=100)
norin.form %>% get_prior(data=norin)
## prior class coef group resp dpar nlpar bound
## (flat) b
## (flat) b MASS
## (flat) b SMR_contr
## (flat) b TRIALHypoxia
## (flat) b TRIALHypoxia:SMR_contr
## (flat) b TRIALLowSalinity
## (flat) b TRIALLowSalinity:SMR_contr
## student_t(3, 16.4, 34.5) Intercept
## student_t(3, 0, 34.5) sd
## student_t(3, 0, 34.5) sd FISHID
## student_t(3, 0, 34.5) sd Intercept FISHID
## student_t(3, 0, 34.5) sigma
## source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## default
norin.brm <- brm(norin.form,
data=norin,
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0)
## Compiling Stan program...
## Start sampling
This tells us:
median(norin$CHANGE)mad(norin$CHANGE)brms uses flat (inproper) priors for all population
effectsmad(norin$CHANGE)The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):
median(norin$CHANGE)mad(norin$CHANGE)sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)sd(norin$CHANGE)sqrt(sd(norin$CHANGE))Note, for hierarchical models, the model will tend to want to have a large \(sigma\) in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:
I will also overlay the raw data for comparison.
norin %>%
group_by(TRIAL) %>%
summarise(median(CHANGE),
mad(CHANGE))
norin %>%
group_by(TRIAL, FISHID) %>%
summarise(median = median(CHANGE),
MAD = mad(CHANGE)) %>%
ungroup(FISHID) %>%
summarise(sd(median))
## `summarise()` has grouped output by 'TRIAL'. You can override using the
## `.groups` argument.
sd(norin$CHANGE)/apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
## TRIALHypoxia TRIALLowSalinity
## 71.597218 71.597218
## SMR_contr MASS
## 54.058070 4.237486
## TRIALHypoxia:SMR_contr TRIALLowSalinity:SMR_contr
## 13.676626 13.676626
standist::visualize("normal(16,35)", xlim=c(-10,100))
standist::visualize("normal(0, 70)", "normal(0, 54)", xlim=c(-200,200))
standist::visualize("gamma(2, 1)", "gamma(35, 1)",
xlim=c(-10,100))
priors <- prior(normal(16, 35), class = 'Intercept') +
prior(normal(0, 70), class = 'b', coef = 'TRIALHypoxia') +
prior(normal(0, 70), class = 'b', coef = 'TRIALLowSalinity') +
prior(normal(0, 54), class = 'b', coef = 'SMR_contr') +
prior(normal(0, 4), class = 'b', coef = 'MASS') +
prior(normal(0, 13), class = 'b') +
prior(gamma(35, 1), class = 'sigma') +
prior(cauchy(0, 2.5), class = 'sd')
norin.form <- bf(CHANGE ~ (1|FISHID) + TRIAL*SMR_contr + MASS,
family = gaussian()
)
norin.brm2 <- brm(norin.form,
data = norin,
prior = priors,
sample_prior = 'only',
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)
## Compiling Stan program...
## Start sampling
Note in the above model, the output may have included a warning message alerting us the presence of divergent transitions. Divergent transitions are an indication that the sampler has encountered poor sampling conditions - the more divergent transitions, the more severe the issue.
Typically, divergent transitions are the result of either:
Accordingly, these divergent transitions can be addressed by either:
norin.brm2 %>%
ggpredict(~SMR_contr*TRIAL) %>%
plot(add.data = TRUE)
## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
norin.brm3 <- update(norin.brm2,
sample_prior = 'yes',
control = list(adapt_delta = 0.99),
refresh = 0)
## The desired updates require recompiling the model
## Compiling Stan program...
## Start sampling
save(norin.brm3, file = '../ws/testing/norin.brm3')
norin.brm3 %>%
ggpredict(~SMR_contr*TRIAL) %>%
plot(add.data = TRUE)
## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
norin.brm3 %>% get_variables()
## [1] "b_Intercept" "b_TRIALHypoxia"
## [3] "b_TRIALLowSalinity" "b_SMR_contr"
## [5] "b_MASS" "b_TRIALHypoxia:SMR_contr"
## [7] "b_TRIALLowSalinity:SMR_contr" "sd_FISHID__Intercept"
## [9] "sigma" "r_FISHID[1,Intercept]"
## [11] "r_FISHID[2,Intercept]" "r_FISHID[3,Intercept]"
## [13] "r_FISHID[4,Intercept]" "r_FISHID[5,Intercept]"
## [15] "r_FISHID[6,Intercept]" "r_FISHID[7,Intercept]"
## [17] "r_FISHID[8,Intercept]" "r_FISHID[9,Intercept]"
## [19] "r_FISHID[10,Intercept]" "r_FISHID[11,Intercept]"
## [21] "r_FISHID[12,Intercept]" "r_FISHID[13,Intercept]"
## [23] "r_FISHID[14,Intercept]" "r_FISHID[15,Intercept]"
## [25] "r_FISHID[16,Intercept]" "r_FISHID[17,Intercept]"
## [27] "r_FISHID[18,Intercept]" "r_FISHID[19,Intercept]"
## [29] "r_FISHID[20,Intercept]" "r_FISHID[21,Intercept]"
## [31] "r_FISHID[22,Intercept]" "r_FISHID[23,Intercept]"
## [33] "r_FISHID[24,Intercept]" "r_FISHID[25,Intercept]"
## [35] "r_FISHID[26,Intercept]" "r_FISHID[27,Intercept]"
## [37] "r_FISHID[28,Intercept]" "r_FISHID[29,Intercept]"
## [39] "r_FISHID[30,Intercept]" "r_FISHID[31,Intercept]"
## [41] "r_FISHID[32,Intercept]" "r_FISHID[33,Intercept]"
## [43] "r_FISHID[34,Intercept]" "r_FISHID[35,Intercept]"
## [45] "r_FISHID[36,Intercept]" "r_FISHID[37,Intercept]"
## [47] "r_FISHID[38,Intercept]" "r_FISHID[39,Intercept]"
## [49] "r_FISHID[40,Intercept]" "r_FISHID[41,Intercept]"
## [51] "r_FISHID[42,Intercept]" "r_FISHID[43,Intercept]"
## [53] "r_FISHID[44,Intercept]" "r_FISHID[45,Intercept]"
## [55] "r_FISHID[46,Intercept]" "r_FISHID[47,Intercept]"
## [57] "r_FISHID[48,Intercept]" "r_FISHID[49,Intercept]"
## [59] "r_FISHID[50,Intercept]" "r_FISHID[51,Intercept]"
## [61] "r_FISHID[52,Intercept]" "r_FISHID[53,Intercept]"
## [63] "r_FISHID[54,Intercept]" "r_FISHID[55,Intercept]"
## [65] "r_FISHID[56,Intercept]" "r_FISHID[57,Intercept]"
## [67] "r_FISHID[58,Intercept]" "r_FISHID[59,Intercept]"
## [69] "r_FISHID[60,Intercept]" "prior_Intercept"
## [71] "prior_b_TRIALHypoxia" "prior_b_TRIALLowSalinity"
## [73] "prior_b_SMR_contr" "prior_b_MASS"
## [75] "prior_b_TRIALHypoxia:SMR_contr" "prior_b_TRIALLowSalinity:SMR_contr"
## [77] "prior_sigma" "prior_sd_FISHID"
## [79] "lp__" "accept_stat__"
## [81] "stepsize__" "treedepth__"
## [83] "n_leapfrog__" "divergent__"
## [85] "energy__"
norin.brm3 %>% hypothesis('TRIALHypoxia=0') %>% plot
norin.brm3 %>% hypothesis('SMR_contr=0') %>% plot
norin.brm3 %>% SUYR_prior_and_posterior()
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
While we are here, we might like to explore a random intercept/slope model
priors <- prior(normal(16, 35), class = 'Intercept') +
prior(normal(0, 70), class = 'b', coef = 'TRIALHypoxia') +
prior(normal(0, 70), class = 'b', coef = 'TRIALLowSalinity') +
prior(normal(0, 54), class = 'b', coef = 'SMR_contr') +
prior(normal(0, 4), class = 'b', coef = 'MASS') +
prior(normal(0, 13), class = 'b') +
prior(gamma(35, 1), class = 'sigma') +
prior(cauchy(0, 2.5), class = 'sd') +
prior(lkj_corr_cholesky(1), class = 'L')
norin.form <- bf(CHANGE ~ (TRIAL|FISHID) + TRIAL*SMR_contr + MASS,
family = gaussian()
)
norin.brm4 <- brm(norin.form,
data = norin,
prior = priors,
sample_prior = 'yes',
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0,
control = list(adapt_delta=0.99)
)
## Compiling Stan program...
## Start sampling
save(norin.brm4, file = '../ws/testing/norin.brm4')
norin.brm4 %>% get_variables()
## [1] "b_Intercept"
## [2] "b_TRIALHypoxia"
## [3] "b_TRIALLowSalinity"
## [4] "b_SMR_contr"
## [5] "b_MASS"
## [6] "b_TRIALHypoxia:SMR_contr"
## [7] "b_TRIALLowSalinity:SMR_contr"
## [8] "sd_FISHID__Intercept"
## [9] "sd_FISHID__TRIALHypoxia"
## [10] "sd_FISHID__TRIALLowSalinity"
## [11] "cor_FISHID__Intercept__TRIALHypoxia"
## [12] "cor_FISHID__Intercept__TRIALLowSalinity"
## [13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
## [14] "sigma"
## [15] "r_FISHID[1,Intercept]"
## [16] "r_FISHID[2,Intercept]"
## [17] "r_FISHID[3,Intercept]"
## [18] "r_FISHID[4,Intercept]"
## [19] "r_FISHID[5,Intercept]"
## [20] "r_FISHID[6,Intercept]"
## [21] "r_FISHID[7,Intercept]"
## [22] "r_FISHID[8,Intercept]"
## [23] "r_FISHID[9,Intercept]"
## [24] "r_FISHID[10,Intercept]"
## [25] "r_FISHID[11,Intercept]"
## [26] "r_FISHID[12,Intercept]"
## [27] "r_FISHID[13,Intercept]"
## [28] "r_FISHID[14,Intercept]"
## [29] "r_FISHID[15,Intercept]"
## [30] "r_FISHID[16,Intercept]"
## [31] "r_FISHID[17,Intercept]"
## [32] "r_FISHID[18,Intercept]"
## [33] "r_FISHID[19,Intercept]"
## [34] "r_FISHID[20,Intercept]"
## [35] "r_FISHID[21,Intercept]"
## [36] "r_FISHID[22,Intercept]"
## [37] "r_FISHID[23,Intercept]"
## [38] "r_FISHID[24,Intercept]"
## [39] "r_FISHID[25,Intercept]"
## [40] "r_FISHID[26,Intercept]"
## [41] "r_FISHID[27,Intercept]"
## [42] "r_FISHID[28,Intercept]"
## [43] "r_FISHID[29,Intercept]"
## [44] "r_FISHID[30,Intercept]"
## [45] "r_FISHID[31,Intercept]"
## [46] "r_FISHID[32,Intercept]"
## [47] "r_FISHID[33,Intercept]"
## [48] "r_FISHID[34,Intercept]"
## [49] "r_FISHID[35,Intercept]"
## [50] "r_FISHID[36,Intercept]"
## [51] "r_FISHID[37,Intercept]"
## [52] "r_FISHID[38,Intercept]"
## [53] "r_FISHID[39,Intercept]"
## [54] "r_FISHID[40,Intercept]"
## [55] "r_FISHID[41,Intercept]"
## [56] "r_FISHID[42,Intercept]"
## [57] "r_FISHID[43,Intercept]"
## [58] "r_FISHID[44,Intercept]"
## [59] "r_FISHID[45,Intercept]"
## [60] "r_FISHID[46,Intercept]"
## [61] "r_FISHID[47,Intercept]"
## [62] "r_FISHID[48,Intercept]"
## [63] "r_FISHID[49,Intercept]"
## [64] "r_FISHID[50,Intercept]"
## [65] "r_FISHID[51,Intercept]"
## [66] "r_FISHID[52,Intercept]"
## [67] "r_FISHID[53,Intercept]"
## [68] "r_FISHID[54,Intercept]"
## [69] "r_FISHID[55,Intercept]"
## [70] "r_FISHID[56,Intercept]"
## [71] "r_FISHID[57,Intercept]"
## [72] "r_FISHID[58,Intercept]"
## [73] "r_FISHID[59,Intercept]"
## [74] "r_FISHID[60,Intercept]"
## [75] "r_FISHID[1,TRIALHypoxia]"
## [76] "r_FISHID[2,TRIALHypoxia]"
## [77] "r_FISHID[3,TRIALHypoxia]"
## [78] "r_FISHID[4,TRIALHypoxia]"
## [79] "r_FISHID[5,TRIALHypoxia]"
## [80] "r_FISHID[6,TRIALHypoxia]"
## [81] "r_FISHID[7,TRIALHypoxia]"
## [82] "r_FISHID[8,TRIALHypoxia]"
## [83] "r_FISHID[9,TRIALHypoxia]"
## [84] "r_FISHID[10,TRIALHypoxia]"
## [85] "r_FISHID[11,TRIALHypoxia]"
## [86] "r_FISHID[12,TRIALHypoxia]"
## [87] "r_FISHID[13,TRIALHypoxia]"
## [88] "r_FISHID[14,TRIALHypoxia]"
## [89] "r_FISHID[15,TRIALHypoxia]"
## [90] "r_FISHID[16,TRIALHypoxia]"
## [91] "r_FISHID[17,TRIALHypoxia]"
## [92] "r_FISHID[18,TRIALHypoxia]"
## [93] "r_FISHID[19,TRIALHypoxia]"
## [94] "r_FISHID[20,TRIALHypoxia]"
## [95] "r_FISHID[21,TRIALHypoxia]"
## [96] "r_FISHID[22,TRIALHypoxia]"
## [97] "r_FISHID[23,TRIALHypoxia]"
## [98] "r_FISHID[24,TRIALHypoxia]"
## [99] "r_FISHID[25,TRIALHypoxia]"
## [100] "r_FISHID[26,TRIALHypoxia]"
## [101] "r_FISHID[27,TRIALHypoxia]"
## [102] "r_FISHID[28,TRIALHypoxia]"
## [103] "r_FISHID[29,TRIALHypoxia]"
## [104] "r_FISHID[30,TRIALHypoxia]"
## [105] "r_FISHID[31,TRIALHypoxia]"
## [106] "r_FISHID[32,TRIALHypoxia]"
## [107] "r_FISHID[33,TRIALHypoxia]"
## [108] "r_FISHID[34,TRIALHypoxia]"
## [109] "r_FISHID[35,TRIALHypoxia]"
## [110] "r_FISHID[36,TRIALHypoxia]"
## [111] "r_FISHID[37,TRIALHypoxia]"
## [112] "r_FISHID[38,TRIALHypoxia]"
## [113] "r_FISHID[39,TRIALHypoxia]"
## [114] "r_FISHID[40,TRIALHypoxia]"
## [115] "r_FISHID[41,TRIALHypoxia]"
## [116] "r_FISHID[42,TRIALHypoxia]"
## [117] "r_FISHID[43,TRIALHypoxia]"
## [118] "r_FISHID[44,TRIALHypoxia]"
## [119] "r_FISHID[45,TRIALHypoxia]"
## [120] "r_FISHID[46,TRIALHypoxia]"
## [121] "r_FISHID[47,TRIALHypoxia]"
## [122] "r_FISHID[48,TRIALHypoxia]"
## [123] "r_FISHID[49,TRIALHypoxia]"
## [124] "r_FISHID[50,TRIALHypoxia]"
## [125] "r_FISHID[51,TRIALHypoxia]"
## [126] "r_FISHID[52,TRIALHypoxia]"
## [127] "r_FISHID[53,TRIALHypoxia]"
## [128] "r_FISHID[54,TRIALHypoxia]"
## [129] "r_FISHID[55,TRIALHypoxia]"
## [130] "r_FISHID[56,TRIALHypoxia]"
## [131] "r_FISHID[57,TRIALHypoxia]"
## [132] "r_FISHID[58,TRIALHypoxia]"
## [133] "r_FISHID[59,TRIALHypoxia]"
## [134] "r_FISHID[60,TRIALHypoxia]"
## [135] "r_FISHID[1,TRIALLowSalinity]"
## [136] "r_FISHID[2,TRIALLowSalinity]"
## [137] "r_FISHID[3,TRIALLowSalinity]"
## [138] "r_FISHID[4,TRIALLowSalinity]"
## [139] "r_FISHID[5,TRIALLowSalinity]"
## [140] "r_FISHID[6,TRIALLowSalinity]"
## [141] "r_FISHID[7,TRIALLowSalinity]"
## [142] "r_FISHID[8,TRIALLowSalinity]"
## [143] "r_FISHID[9,TRIALLowSalinity]"
## [144] "r_FISHID[10,TRIALLowSalinity]"
## [145] "r_FISHID[11,TRIALLowSalinity]"
## [146] "r_FISHID[12,TRIALLowSalinity]"
## [147] "r_FISHID[13,TRIALLowSalinity]"
## [148] "r_FISHID[14,TRIALLowSalinity]"
## [149] "r_FISHID[15,TRIALLowSalinity]"
## [150] "r_FISHID[16,TRIALLowSalinity]"
## [151] "r_FISHID[17,TRIALLowSalinity]"
## [152] "r_FISHID[18,TRIALLowSalinity]"
## [153] "r_FISHID[19,TRIALLowSalinity]"
## [154] "r_FISHID[20,TRIALLowSalinity]"
## [155] "r_FISHID[21,TRIALLowSalinity]"
## [156] "r_FISHID[22,TRIALLowSalinity]"
## [157] "r_FISHID[23,TRIALLowSalinity]"
## [158] "r_FISHID[24,TRIALLowSalinity]"
## [159] "r_FISHID[25,TRIALLowSalinity]"
## [160] "r_FISHID[26,TRIALLowSalinity]"
## [161] "r_FISHID[27,TRIALLowSalinity]"
## [162] "r_FISHID[28,TRIALLowSalinity]"
## [163] "r_FISHID[29,TRIALLowSalinity]"
## [164] "r_FISHID[30,TRIALLowSalinity]"
## [165] "r_FISHID[31,TRIALLowSalinity]"
## [166] "r_FISHID[32,TRIALLowSalinity]"
## [167] "r_FISHID[33,TRIALLowSalinity]"
## [168] "r_FISHID[34,TRIALLowSalinity]"
## [169] "r_FISHID[35,TRIALLowSalinity]"
## [170] "r_FISHID[36,TRIALLowSalinity]"
## [171] "r_FISHID[37,TRIALLowSalinity]"
## [172] "r_FISHID[38,TRIALLowSalinity]"
## [173] "r_FISHID[39,TRIALLowSalinity]"
## [174] "r_FISHID[40,TRIALLowSalinity]"
## [175] "r_FISHID[41,TRIALLowSalinity]"
## [176] "r_FISHID[42,TRIALLowSalinity]"
## [177] "r_FISHID[43,TRIALLowSalinity]"
## [178] "r_FISHID[44,TRIALLowSalinity]"
## [179] "r_FISHID[45,TRIALLowSalinity]"
## [180] "r_FISHID[46,TRIALLowSalinity]"
## [181] "r_FISHID[47,TRIALLowSalinity]"
## [182] "r_FISHID[48,TRIALLowSalinity]"
## [183] "r_FISHID[49,TRIALLowSalinity]"
## [184] "r_FISHID[50,TRIALLowSalinity]"
## [185] "r_FISHID[51,TRIALLowSalinity]"
## [186] "r_FISHID[52,TRIALLowSalinity]"
## [187] "r_FISHID[53,TRIALLowSalinity]"
## [188] "r_FISHID[54,TRIALLowSalinity]"
## [189] "r_FISHID[55,TRIALLowSalinity]"
## [190] "r_FISHID[56,TRIALLowSalinity]"
## [191] "r_FISHID[57,TRIALLowSalinity]"
## [192] "r_FISHID[58,TRIALLowSalinity]"
## [193] "r_FISHID[59,TRIALLowSalinity]"
## [194] "r_FISHID[60,TRIALLowSalinity]"
## [195] "prior_Intercept"
## [196] "prior_b_TRIALHypoxia"
## [197] "prior_b_TRIALLowSalinity"
## [198] "prior_b_SMR_contr"
## [199] "prior_b_MASS"
## [200] "prior_b_TRIALHypoxia:SMR_contr"
## [201] "prior_b_TRIALLowSalinity:SMR_contr"
## [202] "prior_sigma"
## [203] "prior_sd_FISHID"
## [204] "prior_cor_FISHID"
## [205] "lp__"
## [206] "accept_stat__"
## [207] "stepsize__"
## [208] "treedepth__"
## [209] "n_leapfrog__"
## [210] "divergent__"
## [211] "energy__"
norin.brm4 %>% hypothesis('TRIALHypoxia=0') %>% plot
norin.brm4 %>% hypothesis('SMR_contr=0') %>% plot
norin.brm4 %>% SUYR_prior_and_posterior()
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
norin.brm4 %>%
posterior_samples %>%
dplyr::select(-`lp__`) %>%
pivot_longer(everything(), names_to = 'key') %>%
filter(!str_detect(key, '^r')) %>%
mutate(Type = ifelse(str_detect(key, 'prior'), 'Prior', 'Posterior'),
Class = case_when(
str_detect(key, '(^b|^prior).*Intercept$') ~ 'Intercept',
str_detect(key, 'b_TRIAL.*|prior_b_TRIAL.*') & str_detect(key, '.*\\:.*') ~ 'TRIAL',
str_detect(key, 'b_SMR_contr|prior_b_SMR_contr') ~ 'SMR',
str_detect(key, 'b_MASS|prior_b_MASS') ~ 'MASS',
str_detect(key, '.*\\:.*|prior_b_.*\\:.*') ~ 'Interaction',
str_detect(key, 'sd') ~ 'sd',
str_detect(key, '^cor|prior_cor') ~ 'cor',
str_detect(key, 'sigma') ~ 'sigma'),
Par = str_replace(key, 'b_', '')) %>%
ggplot(aes(x = Type, y = value, color = Par)) +
stat_pointinterval(position = position_dodge())+
facet_wrap(~Class, scales = 'free')
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
We can compare the two models using LOO
(l.1 <- norin.brm3 %>% loo())
##
## Computed from 2400 by 180 log-likelihood matrix
##
## Estimate SE
## elpd_loo -833.4 12.8
## p_loo 15.2 2.7
## looic 1666.7 25.5
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 177 98.3% 275
## (0.5, 0.7] (ok) 3 1.7% 467
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
(l.2 <- norin.brm4 %>% loo())
## Warning: Found 35 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
##
## Computed from 2400 by 180 log-likelihood matrix
##
## Estimate SE
## elpd_loo -807.0 9.5
## p_loo 54.5 6.1
## looic 1613.9 19.0
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 96 53.3% 522
## (0.5, 0.7] (ok) 49 27.2% 128
## (0.7, 1] (bad) 30 16.7% 21
## (1, Inf) (very bad) 5 2.8% 7
## See help('pareto-k-diagnostic') for details.
loo_compare(l.1, l.2)
## elpd_diff se_diff
## . 0.0 0.0
## . -26.4 6.6
There is substantially more support for the more complex random intercept/slope model over the simpler random intercept only model. Consequently, we will continue with the random intercept/slope model.
In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.
There are a wide variety of tests that range from the big picture, overal chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.
The bayesplot package offers a range of MCMC diagnostics
as well as Posterior Probability Checks (PPC), all of which have a
convenient plot() interface. Lets start with the MCMC
diagnostics.
available_mcmc()
## bayesplot MCMC module:
## mcmc_acf
## mcmc_acf_bar
## mcmc_areas
## mcmc_areas_data
## mcmc_areas_ridges
## mcmc_areas_ridges_data
## mcmc_combo
## mcmc_dens
## mcmc_dens_chains
## mcmc_dens_chains_data
## mcmc_dens_overlay
## mcmc_hex
## mcmc_hist
## mcmc_hist_by_chain
## mcmc_intervals
## mcmc_intervals_data
## mcmc_neff
## mcmc_neff_data
## mcmc_neff_hist
## mcmc_nuts_acceptance
## mcmc_nuts_divergence
## mcmc_nuts_energy
## mcmc_nuts_stepsize
## mcmc_nuts_treedepth
## mcmc_pairs
## mcmc_parcoord
## mcmc_parcoord_data
## mcmc_rank_hist
## mcmc_rank_overlay
## mcmc_recover_hist
## mcmc_recover_intervals
## mcmc_recover_scatter
## mcmc_rhat
## mcmc_rhat_data
## mcmc_rhat_hist
## mcmc_scatter
## mcmc_trace
## mcmc_trace_data
## mcmc_trace_highlight
## mcmc_violin
Of these, we will focus on:
norin.brm4 %>% mcmc_plot(type='trace')
## No divergences to plot.
The chains appear well mixed and very similar
norin.brm4 %>% mcmc_plot(type='acf_bar')
There IS evidence of autocorrelation in the MCMC samples. We should consider thinning further!
norin.brm4 %>% mcmc_plot(type='rhat_hist')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
All Rhat values are below 1.05, suggesting the chains have converged.
neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
norin.brm4 %>% mcmc_plot(type='neff_hist')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
There are a small number of parameters that have very low effective sample sizes. It might be worth exploring the cause(s) of this to determine whether it is concerning.
norin.brm4 %>% mcmc_plot(type='combo')
norin.brm4 %>% mcmc_plot(type='violin')
The rstan package offers a range of MCMC diagnostics.
Lets start with the MCMC diagnostics.
Of these, we will focus on:
norin.brm4 %>% get_variables()
## [1] "b_Intercept"
## [2] "b_TRIALHypoxia"
## [3] "b_TRIALLowSalinity"
## [4] "b_SMR_contr"
## [5] "b_MASS"
## [6] "b_TRIALHypoxia:SMR_contr"
## [7] "b_TRIALLowSalinity:SMR_contr"
## [8] "sd_FISHID__Intercept"
## [9] "sd_FISHID__TRIALHypoxia"
## [10] "sd_FISHID__TRIALLowSalinity"
## [11] "cor_FISHID__Intercept__TRIALHypoxia"
## [12] "cor_FISHID__Intercept__TRIALLowSalinity"
## [13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
## [14] "sigma"
## [15] "r_FISHID[1,Intercept]"
## [16] "r_FISHID[2,Intercept]"
## [17] "r_FISHID[3,Intercept]"
## [18] "r_FISHID[4,Intercept]"
## [19] "r_FISHID[5,Intercept]"
## [20] "r_FISHID[6,Intercept]"
## [21] "r_FISHID[7,Intercept]"
## [22] "r_FISHID[8,Intercept]"
## [23] "r_FISHID[9,Intercept]"
## [24] "r_FISHID[10,Intercept]"
## [25] "r_FISHID[11,Intercept]"
## [26] "r_FISHID[12,Intercept]"
## [27] "r_FISHID[13,Intercept]"
## [28] "r_FISHID[14,Intercept]"
## [29] "r_FISHID[15,Intercept]"
## [30] "r_FISHID[16,Intercept]"
## [31] "r_FISHID[17,Intercept]"
## [32] "r_FISHID[18,Intercept]"
## [33] "r_FISHID[19,Intercept]"
## [34] "r_FISHID[20,Intercept]"
## [35] "r_FISHID[21,Intercept]"
## [36] "r_FISHID[22,Intercept]"
## [37] "r_FISHID[23,Intercept]"
## [38] "r_FISHID[24,Intercept]"
## [39] "r_FISHID[25,Intercept]"
## [40] "r_FISHID[26,Intercept]"
## [41] "r_FISHID[27,Intercept]"
## [42] "r_FISHID[28,Intercept]"
## [43] "r_FISHID[29,Intercept]"
## [44] "r_FISHID[30,Intercept]"
## [45] "r_FISHID[31,Intercept]"
## [46] "r_FISHID[32,Intercept]"
## [47] "r_FISHID[33,Intercept]"
## [48] "r_FISHID[34,Intercept]"
## [49] "r_FISHID[35,Intercept]"
## [50] "r_FISHID[36,Intercept]"
## [51] "r_FISHID[37,Intercept]"
## [52] "r_FISHID[38,Intercept]"
## [53] "r_FISHID[39,Intercept]"
## [54] "r_FISHID[40,Intercept]"
## [55] "r_FISHID[41,Intercept]"
## [56] "r_FISHID[42,Intercept]"
## [57] "r_FISHID[43,Intercept]"
## [58] "r_FISHID[44,Intercept]"
## [59] "r_FISHID[45,Intercept]"
## [60] "r_FISHID[46,Intercept]"
## [61] "r_FISHID[47,Intercept]"
## [62] "r_FISHID[48,Intercept]"
## [63] "r_FISHID[49,Intercept]"
## [64] "r_FISHID[50,Intercept]"
## [65] "r_FISHID[51,Intercept]"
## [66] "r_FISHID[52,Intercept]"
## [67] "r_FISHID[53,Intercept]"
## [68] "r_FISHID[54,Intercept]"
## [69] "r_FISHID[55,Intercept]"
## [70] "r_FISHID[56,Intercept]"
## [71] "r_FISHID[57,Intercept]"
## [72] "r_FISHID[58,Intercept]"
## [73] "r_FISHID[59,Intercept]"
## [74] "r_FISHID[60,Intercept]"
## [75] "r_FISHID[1,TRIALHypoxia]"
## [76] "r_FISHID[2,TRIALHypoxia]"
## [77] "r_FISHID[3,TRIALHypoxia]"
## [78] "r_FISHID[4,TRIALHypoxia]"
## [79] "r_FISHID[5,TRIALHypoxia]"
## [80] "r_FISHID[6,TRIALHypoxia]"
## [81] "r_FISHID[7,TRIALHypoxia]"
## [82] "r_FISHID[8,TRIALHypoxia]"
## [83] "r_FISHID[9,TRIALHypoxia]"
## [84] "r_FISHID[10,TRIALHypoxia]"
## [85] "r_FISHID[11,TRIALHypoxia]"
## [86] "r_FISHID[12,TRIALHypoxia]"
## [87] "r_FISHID[13,TRIALHypoxia]"
## [88] "r_FISHID[14,TRIALHypoxia]"
## [89] "r_FISHID[15,TRIALHypoxia]"
## [90] "r_FISHID[16,TRIALHypoxia]"
## [91] "r_FISHID[17,TRIALHypoxia]"
## [92] "r_FISHID[18,TRIALHypoxia]"
## [93] "r_FISHID[19,TRIALHypoxia]"
## [94] "r_FISHID[20,TRIALHypoxia]"
## [95] "r_FISHID[21,TRIALHypoxia]"
## [96] "r_FISHID[22,TRIALHypoxia]"
## [97] "r_FISHID[23,TRIALHypoxia]"
## [98] "r_FISHID[24,TRIALHypoxia]"
## [99] "r_FISHID[25,TRIALHypoxia]"
## [100] "r_FISHID[26,TRIALHypoxia]"
## [101] "r_FISHID[27,TRIALHypoxia]"
## [102] "r_FISHID[28,TRIALHypoxia]"
## [103] "r_FISHID[29,TRIALHypoxia]"
## [104] "r_FISHID[30,TRIALHypoxia]"
## [105] "r_FISHID[31,TRIALHypoxia]"
## [106] "r_FISHID[32,TRIALHypoxia]"
## [107] "r_FISHID[33,TRIALHypoxia]"
## [108] "r_FISHID[34,TRIALHypoxia]"
## [109] "r_FISHID[35,TRIALHypoxia]"
## [110] "r_FISHID[36,TRIALHypoxia]"
## [111] "r_FISHID[37,TRIALHypoxia]"
## [112] "r_FISHID[38,TRIALHypoxia]"
## [113] "r_FISHID[39,TRIALHypoxia]"
## [114] "r_FISHID[40,TRIALHypoxia]"
## [115] "r_FISHID[41,TRIALHypoxia]"
## [116] "r_FISHID[42,TRIALHypoxia]"
## [117] "r_FISHID[43,TRIALHypoxia]"
## [118] "r_FISHID[44,TRIALHypoxia]"
## [119] "r_FISHID[45,TRIALHypoxia]"
## [120] "r_FISHID[46,TRIALHypoxia]"
## [121] "r_FISHID[47,TRIALHypoxia]"
## [122] "r_FISHID[48,TRIALHypoxia]"
## [123] "r_FISHID[49,TRIALHypoxia]"
## [124] "r_FISHID[50,TRIALHypoxia]"
## [125] "r_FISHID[51,TRIALHypoxia]"
## [126] "r_FISHID[52,TRIALHypoxia]"
## [127] "r_FISHID[53,TRIALHypoxia]"
## [128] "r_FISHID[54,TRIALHypoxia]"
## [129] "r_FISHID[55,TRIALHypoxia]"
## [130] "r_FISHID[56,TRIALHypoxia]"
## [131] "r_FISHID[57,TRIALHypoxia]"
## [132] "r_FISHID[58,TRIALHypoxia]"
## [133] "r_FISHID[59,TRIALHypoxia]"
## [134] "r_FISHID[60,TRIALHypoxia]"
## [135] "r_FISHID[1,TRIALLowSalinity]"
## [136] "r_FISHID[2,TRIALLowSalinity]"
## [137] "r_FISHID[3,TRIALLowSalinity]"
## [138] "r_FISHID[4,TRIALLowSalinity]"
## [139] "r_FISHID[5,TRIALLowSalinity]"
## [140] "r_FISHID[6,TRIALLowSalinity]"
## [141] "r_FISHID[7,TRIALLowSalinity]"
## [142] "r_FISHID[8,TRIALLowSalinity]"
## [143] "r_FISHID[9,TRIALLowSalinity]"
## [144] "r_FISHID[10,TRIALLowSalinity]"
## [145] "r_FISHID[11,TRIALLowSalinity]"
## [146] "r_FISHID[12,TRIALLowSalinity]"
## [147] "r_FISHID[13,TRIALLowSalinity]"
## [148] "r_FISHID[14,TRIALLowSalinity]"
## [149] "r_FISHID[15,TRIALLowSalinity]"
## [150] "r_FISHID[16,TRIALLowSalinity]"
## [151] "r_FISHID[17,TRIALLowSalinity]"
## [152] "r_FISHID[18,TRIALLowSalinity]"
## [153] "r_FISHID[19,TRIALLowSalinity]"
## [154] "r_FISHID[20,TRIALLowSalinity]"
## [155] "r_FISHID[21,TRIALLowSalinity]"
## [156] "r_FISHID[22,TRIALLowSalinity]"
## [157] "r_FISHID[23,TRIALLowSalinity]"
## [158] "r_FISHID[24,TRIALLowSalinity]"
## [159] "r_FISHID[25,TRIALLowSalinity]"
## [160] "r_FISHID[26,TRIALLowSalinity]"
## [161] "r_FISHID[27,TRIALLowSalinity]"
## [162] "r_FISHID[28,TRIALLowSalinity]"
## [163] "r_FISHID[29,TRIALLowSalinity]"
## [164] "r_FISHID[30,TRIALLowSalinity]"
## [165] "r_FISHID[31,TRIALLowSalinity]"
## [166] "r_FISHID[32,TRIALLowSalinity]"
## [167] "r_FISHID[33,TRIALLowSalinity]"
## [168] "r_FISHID[34,TRIALLowSalinity]"
## [169] "r_FISHID[35,TRIALLowSalinity]"
## [170] "r_FISHID[36,TRIALLowSalinity]"
## [171] "r_FISHID[37,TRIALLowSalinity]"
## [172] "r_FISHID[38,TRIALLowSalinity]"
## [173] "r_FISHID[39,TRIALLowSalinity]"
## [174] "r_FISHID[40,TRIALLowSalinity]"
## [175] "r_FISHID[41,TRIALLowSalinity]"
## [176] "r_FISHID[42,TRIALLowSalinity]"
## [177] "r_FISHID[43,TRIALLowSalinity]"
## [178] "r_FISHID[44,TRIALLowSalinity]"
## [179] "r_FISHID[45,TRIALLowSalinity]"
## [180] "r_FISHID[46,TRIALLowSalinity]"
## [181] "r_FISHID[47,TRIALLowSalinity]"
## [182] "r_FISHID[48,TRIALLowSalinity]"
## [183] "r_FISHID[49,TRIALLowSalinity]"
## [184] "r_FISHID[50,TRIALLowSalinity]"
## [185] "r_FISHID[51,TRIALLowSalinity]"
## [186] "r_FISHID[52,TRIALLowSalinity]"
## [187] "r_FISHID[53,TRIALLowSalinity]"
## [188] "r_FISHID[54,TRIALLowSalinity]"
## [189] "r_FISHID[55,TRIALLowSalinity]"
## [190] "r_FISHID[56,TRIALLowSalinity]"
## [191] "r_FISHID[57,TRIALLowSalinity]"
## [192] "r_FISHID[58,TRIALLowSalinity]"
## [193] "r_FISHID[59,TRIALLowSalinity]"
## [194] "r_FISHID[60,TRIALLowSalinity]"
## [195] "prior_Intercept"
## [196] "prior_b_TRIALHypoxia"
## [197] "prior_b_TRIALLowSalinity"
## [198] "prior_b_SMR_contr"
## [199] "prior_b_MASS"
## [200] "prior_b_TRIALHypoxia:SMR_contr"
## [201] "prior_b_TRIALLowSalinity:SMR_contr"
## [202] "prior_sigma"
## [203] "prior_sd_FISHID"
## [204] "prior_cor_FISHID"
## [205] "lp__"
## [206] "accept_stat__"
## [207] "stepsize__"
## [208] "treedepth__"
## [209] "n_leapfrog__"
## [210] "divergent__"
## [211] "energy__"
pars <- norin.brm4 %>% get_variables()
pars <- str_extract(pars, '^b_.*|^sigma$|^sd.*') %>% na.omit()
norin.brm4$fit %>%
stan_trace(pars = pars)
The chains appear well mixed and very similar
norin.brm4$fit %>%
stan_ac(pars = pars)
There IS evidence of autocorrelation in the MCMC samples. We should consider thinning further!
norin.brm4$fit %>% stan_rhat()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
norin.brm4$fit %>% stan_ess()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
There are a small number of parameters that have very low effective sample sizes. It might be worth exploring the cause(s) of this to determine whether it is concerning.
norin.brm4$fit %>%
stan_dens(separate_chains = TRUE, pars = pars)
The ggmean package also as a set of MCMC diagnostic
functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
norin.ggs <- norin.brm4 %>% ggs(burnin = FALSE, inc_warmup = FALSE)
## Warning in custom.sort(D$Parameter): NAs introduced by coercion
norin.ggs %>% ggs_traceplot()
The chains appear well mixed and very similar
ggs_autocorrelation(norin.ggs)
There is no evidence of autocorrelation in the MCMC samples
ggs_Rhat(norin.ggs, scaling = 1.01)
All Rhat values are below 1.01, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
ggs_effective(norin.ggs)
Ratios all very high.
ggs_crosscorrelation(norin.ggs)
ggs_grb(norin.ggs)
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
available_ppc()
## bayesplot PPC module:
## ppc_bars
## ppc_bars_grouped
## ppc_boxplot
## ppc_data
## ppc_dens
## ppc_dens_overlay
## ppc_dens_overlay_grouped
## ppc_ecdf_overlay
## ppc_ecdf_overlay_grouped
## ppc_error_binned
## ppc_error_hist
## ppc_error_hist_grouped
## ppc_error_scatter
## ppc_error_scatter_avg
## ppc_error_scatter_avg_vs_x
## ppc_freqpoly
## ppc_freqpoly_grouped
## ppc_hist
## ppc_intervals
## ppc_intervals_data
## ppc_intervals_grouped
## ppc_km_overlay
## ppc_loo_intervals
## ppc_loo_pit
## ppc_loo_pit_data
## ppc_loo_pit_overlay
## ppc_loo_pit_qq
## ppc_loo_ribbon
## ppc_ribbon
## ppc_ribbon_data
## ppc_ribbon_grouped
## ppc_rootogram
## ppc_scatter
## ppc_scatter_avg
## ppc_scatter_avg_grouped
## ppc_stat
## ppc_stat_2d
## ppc_stat_freqpoly_grouped
## ppc_stat_grouped
## ppc_violin_grouped
norin.brm4 %>% pp_check(type = 'dens_overlay', ndraws = 100)
The model draws appear to be consistent with the observed data.
norin.brm4 %>% pp_check(type = 'error_scatter_avg')
## Using all posterior draws for ppc type 'error_scatter_avg' by default.
This is not really interpretable
norin.brm4 %>% pp_check(group = 'TRIAL', type = 'intervals')
## Using all posterior draws for ppc type 'intervals' by default.
## Warning: The following arguments were unrecognized and ignored: group
norin.brm3 %>% pp_check(group = 'TRIAL', type = 'intervals_grouped')
## Using all posterior draws for ppc type 'intervals_grouped' by default.
norin.brm3 %>% pp_check(group = 'TRIAL', type = 'violin_grouped')
## Using all posterior draws for ppc type 'violin_grouped' by default.
The shinystan package allows the full suite of MCMC
diagnostics and posterior predictive checks to be accessed via a web
interface.
#library(shinystan)
#launch_shinystan(norin.brm2)
DHARMa residuals provide very useful diagnostics. Unfortunately, we
cannot directly use the simulateResiduals() function to
generate the simulated residuals. However, if we are willing to
calculate some of the components ourself, we can still obtain the
simulated residuals from the fitted stan model.
We need to supply:
preds <- norin.brm4 %>% posterior_predict(ndraws = 250, summary = FALSE)
norin.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = norin$CHANGE,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = FALSE)
plot(norin.resids, quantreg = FALSE)
Conclusions:
norin.brm4 %>%
conditional_effects() %>%
plot(ask = FALSE, points = TRUE, plot = FALSE) %>%
wrap_plots()
norin.brm4 %>%
ggpredict(terms = c("SMR_contr", "TRIAL")) %>%
plot(add.data = TRUE)
## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
norin.brm4 %>%
ggemmeans(~SMR_contr*TRIAL) %>%
plot(add.data = TRUE)
It is not really possible to do this via the fitted draws as it would
not be marginalising over MASS or FISHID.
brms captures the MCMC samples from stan
within the returned list. There are numerous ways to retrieve and
summarise these samples. The first three provide convenient numeric
summaries from which you can draw conclusions, the last four provide
ways of obtaining the full posteriors.
The summary() method generates simple summaries (mean,
standard deviation as well as 10, 50 and 90 percentiles).
norin.brm4 %>% summary()
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + MASS
## Data: norin (Number of observations: 180)
## Draws: 3 chains, each with iter = 1000; warmup = 200; thin = 5;
## total post-warmup draws = 480
##
## Group-Level Effects:
## ~FISHID (Number of levels: 60)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 5.11 3.48 0.19 12.31 1.00
## sd(TRIALHypoxia) 2.46 2.22 0.08 8.34 1.00
## sd(TRIALLowSalinity) 26.24 4.50 17.65 35.41 1.00
## cor(Intercept,TRIALHypoxia) -0.11 0.50 -0.91 0.85 1.00
## cor(Intercept,TRIALLowSalinity) 0.27 0.39 -0.55 0.92 1.00
## cor(TRIALHypoxia,TRIALLowSalinity) -0.19 0.47 -0.93 0.76 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 1086 1792
## sd(TRIALHypoxia) 1831 2066
## sd(TRIALLowSalinity) 1577 1985
## cor(Intercept,TRIALHypoxia) 2120 2210
## cor(Intercept,TRIALLowSalinity) 727 812
## cor(TRIALHypoxia,TRIALLowSalinity) 623 1582
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 140.76 21.78 98.54 184.66 1.00 2262
## TRIALHypoxia -91.92 24.06 -136.82 -44.63 1.00 2368
## TRIALLowSalinity -60.10 30.53 -119.28 1.01 1.00 2079
## SMR_contr -18.13 3.65 -25.37 -10.87 1.00 2259
## MASS 0.09 0.23 -0.37 0.55 1.00 2419
## TRIALHypoxia:SMR_contr 7.99 4.60 -1.25 16.53 1.00 2287
## TRIALLowSalinity:SMR_contr 3.55 5.85 -8.03 15.03 1.00 2024
## Tail_ESS
## Intercept 2452
## TRIALHypoxia 2079
## TRIALLowSalinity 2326
## SMR_contr 2283
## MASS 2369
## TRIALHypoxia:SMR_contr 2325
## TRIALLowSalinity:SMR_contr 2290
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 18.14 1.73 14.62 21.45 1.00 1380 1553
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Conclusions:
the intercept indicates that the average metabolic change associated with the High Temperature trial (at the mean SMR_contr and MASS) across all fish was 140.76.
for the same SMR_contr and MASS, the average metabolic change was 91.92 less in the Hypoxia trial and 60.1 less in the Salinity trial than the High Temperature trial.
in the High Temperature trial, every 1 unit increase in SMR_contr
was associated with a decline in metabolic change of
round(-1*norin.sum$fixed[4,1], 2) units.
in the High Temperature trial (and mean SMR_contr), every 1 unit
increase in MASS was associated with an increase in metabolic change of
round(-1*norin.sum$fixed[5,1], 2) units.
the rate of change in metabolic change associated with the
Hypoxia trial was round(-1*norin.sum$fixed[6,1], 2) units
higher than that of the High Temperature trial and
round(-1*norin.sum$fixed[7,1], 2) units higher in the Low
Salinity trial.
the variance in intercepts across all fish is 5.11
the variance in responses between High Temperature and Hypoxia trials across all fish is 2.46.
the variance in responses between High Temperature and Low Salinity trials across all fish is 26.24.
the scale of variance between fish within a treatment
(sigma, 18.14) is substantially higher than that both within the High
Temperature trial and between the High Temperature and Hypoxia trials,
yet lower than that between High Temperature and Low Salinity
trials.
norin.brm4 %>% as_draws_df()